home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / BESSJ.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  1KB  |  52 lines
  1. FUNCTION bessj(n: integer; x: real): real;
  2. CONST
  3.    iacc=40;
  4.    bigno=1.0e10;
  5.    bigni=1.0e-10;
  6. VAR
  7.    bj,bjm,bjp,sum,tox,ans: real;
  8.    j,jsum,m: integer;
  9. BEGIN
  10.    IF (n < 2) THEN BEGIN
  11.       writeln('pause in BESSJ'); readln
  12.    END;
  13.    IF (x=0.0) THEN ans := 0.0
  14.    ELSE IF (abs(x) > 1.0*n) THEN BEGIN
  15.       tox := 2.0/abs(x);
  16.       bjm := bessj0(abs(x));
  17.       bj := bessj1(abs(x));
  18.       FOR j := 1 TO n-1 DO BEGIN
  19.          bjp := j*tox*bj-bjm;
  20.          bjm := bj;
  21.          bj := bjp
  22.       END;
  23.       ans := bj
  24.    END ELSE BEGIN
  25.       tox := 2.0/abs(x);
  26.       m := 2*((n+trunc(sqrt(1.0*(iacc*n)))) DIV 2);
  27.       ans := 0.0;
  28.       jsum := 0;
  29.       sum := 0.0;
  30.       bjp := 0.0;
  31.       bj := 1.0;
  32.       FOR j := m DOWNTO 1 DO BEGIN
  33.          bjm := j*tox*bj-bjp;
  34.          bjp := bj;
  35.          bj := bjm;
  36.          IF (abs(bj) > bigno) THEN BEGIN
  37.             bj := bj*bigni;
  38.             bjp := bjp*bigni;
  39.             ans := ans*bigni;
  40.             sum := sum*bigni
  41.          END;
  42.          IF (jsum <> 0) THEN sum := sum+bj;
  43.          jsum := 1-jsum;
  44.          IF (j = n) THEN ans := bjp
  45.       END;
  46.       sum := 2.0*sum-bj;
  47.       ans := ans/sum
  48.    END;
  49.    IF (x<0.0) AND ((n MOD 2)=1) THEN ans := -ans;
  50.    bessj := ans
  51. END;
  52.